arXiv: 1507.0077 lv2 [astro-ph.GA] 29 Nov 2015 


Mon. Not. R. Astron. Soc. 000, 000-000 (0000) Printed 1 December 2015 (MN style file v2.2) 


The orbital PDF: the dynamical state of Milky Way sized 
haloes and the intrinsic uncertainty in the determination of 
their masses 

Jiaxin Han, 1 *, Wenting Wang, 1 , Shaun Cole, 1 , Carlos S. Frenk, 1 

1 Institute for Computational Cosmology, Department of Physics, University of Durham, South Road, Durham, DH1 3LE 


1 December 2015 


ABSTRACT 

Using realistic cosmological simulations of Milky Way sized haloes, we study their 
dynamical state and the accuracy of inferring their mass profiles with steady-state 
models of dynamical tracers. We use a new method that describes the phase-space 
distribution of a steady-state tracer population in a spherical potential without any 
assumption regarding the distribution of their orbits. Applying the method to five 
haloes from the Aquarius ACDM N-body simulation, we find that dark matter parti¬ 
cles are an accurate tracer that enables the halo mass and concentration parameters 
to be recovered with an accuracy of 5%. Assuming a potential profile of the NFW form 
does not significantly affect the fits in most cases, except for halo A whose density 
profile differs significantly from the NFW form, leading to a 30% bias in the dynam¬ 
ically fitted parameters. The existence of substructures in the dark matter tracers 
only affects the fits by ^ 1%. Applying the method to mock stellar haloes generated 
by a particle-tagging technique, we find the stars are farther from equilibrium than 
dark matter particles, yielding a systematic bias of ^ 20% in the inferred mass and 
concentration parameter. The level of systematic biases obtained from a conventional 
distribution function fit to stars is comparable to ours, while similar fits to DM tracers 
are significantly biased in contrast to our fits. In line with previous studies, the mass 
bias is much reduced near the tracer half-mass radius. 

Key words: dark matter - galaxies: haloes - galaxies: kinematics and dynamics - 
Galaxy: fundamental parameters - methods: data analysis 


1 INTRODUCTION 

Dynamical modelling is of fundamental importance in the 
determination of the mass distribution of dark matter 
haloes. To constrain the total mass distribution or the grav¬ 
itational potential, a large family of dynamical methods 
work by fitting a proposed potential-dependent distribution 
function (DF) to the observed phase-space distribution of 
a tracer population. In such modelling, one should make as 
few assumptions as possible so as to avoid biasing the re¬ 
sults. In practice, a required minimal assumption is that the 
system is in steady state, so that modelling the tracer DF 
with a single observational snapshot is informative without 
requiring the observation to take place at any special mo¬ 
ment. However, most existing methods involve additional 
assumptions, for example, about the distribution of orbits, 
the functional form of the DF, or the spatial distribution 
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of tracer partic les outside the o bservational window. In a 
previous paper (|Han et al.ll2015l . hereafter Paper I), we de¬ 
veloped a method that can be used to infer the potential 
while only making the assumption that the tracer popula¬ 
tion is in a steady state. In particular, taking a spherical po¬ 
tential as an example, we have shown that the steady-state 
property translates into a fundamental orbital Probability 
Density Function (oPDF), which provides enough informa¬ 
tion to enable the inference of the halo potential. Applying 
this method to a set of steady-state tracers in an NFW po¬ 
tential generated from Monte-Carlo simulations, we showed 
that the method is able to recover the true potential. While 
spherical symmetry is assumed, all the steps of the method 
can be generalized to non-spherical cases. 

A realistic halo from cosmological simulations or one in 
the real universe may violate the assumptions of our method 
in several ways. For example, spherical symm etry is only 
approximate since we k now haloes are triaxial (|Frenk et al.1 
Il98fil: lj ing &; Sut3l2002h . Also, the potential and the distri- 
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bution of any tracer are not strictly static as haloes evolve 
with time. Finally, real haloes are not smooth structures, 
since they contain many subhaloes. In this work, we ap¬ 
ply the oPDF method to the dynamical distribution of dark 
matter and star particles in simulated haloes, to explore the 
extent to which the tracers in a real halo satisfy our model 
assumptions. 

One important motivation for this work is to provide 
a generic assessment of what to expect for the accuracy of 
dynamical mass estimates of the Milky Way (MW) halo. 
The mass of the MW plays a crucial role in i nterpret ¬ 
ing many of the Local Group observations (Wang et al.l 
120121 : iKennedv et al.l [20141 : ICautun et al.l [2014 b However, 
dynamically inferred masses in the literature vary widely, 
ranging from 0. 5 x 10 12 to 2.5 x 10 12 M^ across differ- 
ent studies (e.g., Wilkinson & Evans 1999; Xue_e^jl.]2008; 


Gnedin et al 1l20ld: iGibbons et al 1 120141 : 1 Williams fc Evansl 

2015al : see I Wang et al.l 120151 for a recent compilation of 


measurements). At least part of the discrepancy originates 
from the different assumptions involved in different methods. 
Hence, it is interesting to investigate the intrinsic accuracy 
of a generic dynamical method that makes minimal assump¬ 
tions, which could then be interpreted as a lower limit on 
the systematic uncertainty in dynamical mass estimation. 
Such a study is also timely given the huge amount of phase- 
space data for stars in the Galaxy being ob tained b y a new 
gener ation of instruments such as GAIA (| Perryman et al.l 

l200ih . 


To this end, we apply our generic dynamical method to 
five haloes from the Aquarius simulations, a set of cosmolog¬ 
ical zoom-in simulations of the formation a nd evolution of 
MW sized haloes in the ACDM cosmology (|Springel et al.l 
2008). We fit for the mass and concentration parameters of 
each halo using both the dark matter particles, and the “halo 
stars” from the particle tagging method of ICooper et al.l 
(1201(f ) as tracers. We find that while the dark matter (DM) 
tracers recover the halo parameters accurately, the tagged 
stars result in ~ 20% bias in the dynamically fitted parame¬ 
ters. We give a brief review of the oPDF method in Section[2] 
The applications to DM and stars are presented in Sections^ 
and[5j with the data described in Section [3] and a discussion 
on the half-mass constraint in Section [6] We summarize the 
results and conclude in Section 0 


2 THE OPDF METHOD 

Below we briefly review the oPDF method developed in 
Paper I. A likelihood estimator and a non-parametric pro¬ 
file reconstruction method were developed in Paper I which 
show similar efficiency in making use of the dynamical infor¬ 
mation. We restrict our attention to the likelihood method 
throughout this paper. 


2.1 The oPDF 

In a steady-state system, phase space continuity implies a 
fundamental DF, 

dP(A| orbit)/dA oc d£(A|orbit)/dA, (1) 

where A is an affine parameter specifying the position of a 
particle on a given orbit. That is, for any given orbit, the 


probability of observing a particle at a given position A is 
proportional to the time it spends at that position. In a 
spherical potential, the orbits of particles are described by 
their conserved binding energy, E — — + v /) + ip(r)) 

and conserved angular momentum, L = rv t , where v T and 
vt are the radial and tangential velocities, and ip(r) is the 
potential at radius r. Taking r as the affine parameter, Equa¬ 
tion © becomes, 

dt 1 dr 

dp(rlE ' L) = 7di = TW\’ (2) 


where T = dr/\v Y \ is the period of half an orbit, with 
r p and r a being the peri- and apo-centre radii of the orbit. 
When radial cuts (r m in,r max ) are imposed, we only need 
to replace the orbital limits, r a , with min(r a , r ma x) and r p 
with max(r p , r m in), since Equation [I] holds within any radial 
range. Taking the radial action angle, 0 , which we call the 
phase angle, as the affine parameter, the oPDF becomes a 
uniform distribution, 


where 


dP(0\Ej L) = d<9, 



( 3 ) 

( 4 ) 


This uniform distribution with 6 £ [0,1] is also 

known as the random phase principle or orbital 
roulette ([Beloborodov fc Levinlfiool b 


2.2 Uniform phase diagnostics 

For a steady sta te tracer, if one defines a no rmalized mean 
phase deviation ([Beloborodov &; Levinll2004h by 

0 = v / T2iV(<9-0.5), (5) 

then when the sample size, N, is large enough the uniform 
phase distribution of 6 should result in © being distributed 
like a standard normal variable. Hence, for a real sample, 
@ 2 can be used as a measure of the difference of the actual 
phase distribution from the expected uniform distribution. 


2.3 The radial likelihood estimator 

Given a tracer and an assumed potential, one can predict 
the expected radial PDF of each tracer particle using 

1 N 

dP(r) = -J2 dP ( r ( 6 ) 

3 = 1 


where Ej and Lj are the energy and angular momentum of 
particle j under the assumed potential. If we bin the data 
radially into m bins, the expected number of particles in the 
z-th bin is given by 


fa = N 



d P{r 
dr 


-dr, 


( 7 ) 
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where r\^ and r Uj i are the lower and upper bin edges. The 
binned radial likelihood is given by: 


m 


n«r ex p( _ ”i) 

i=1 

(8) 

m 

exp(— jv) H^r% 

(9) 


i= 1 


where rti is the observed number of particles in the z-th 
bin. The best-fitting potential is defined to be the one that 
maximizes this likelihood. 


3 DATA 


We use the Aquarius simulations ([Springel et al.l [2QQ8h . a 
set of cosmological zoom-in simulations of the formation and 
evolution of MW-sized haloes, for this analysis. The five sim¬ 
ulations we use (labelled “A” to “E”) were run at a series of 
resolutions and we only use the second highest resolution 
(level- 2 ) runs, which have a particle mass of ~ 1 O 4 M 0 so 
that each halo is resolved with ~ 10 8 particles. We consider 
two types of tracers in the halo: DM particles and star par¬ 
ticles. Because Aquarius is a DM-only simulation, the star 
particles are a subs et of DM particles se lected with a particle 
tagging technique ([Cooper et al.l[ 2010 ). 

The oPDF method laid out above assumes a steady- 
state system with a spherical potential. The real halo may 
deviate from these assumptions in many respects, for exam¬ 
ple, by being aspherical, evolving or having substructures. 
We expect these deviations to bias the fit, and our aim is 
to quantify these systematic errors. To this end, we will use 
a large sample to ensure that the statistical noise as in¬ 
ferred from the likelihood estimator is much smaller than 
the level of accuracy of interest. Using Monte Carlo realiza¬ 
tions we found in Paper I that the typical error in halo profile 
parameters is 0.1/ yj iV/1000 dex for N particles. Wherever 
possible, we will use samples with N ~ 10 6 particles lead¬ 
ing to statistical errors of the order of only ~ 1 percent in 
the dynamically derived parameters, the mass, M, and the 
concentration, c. 


3.1 DM Samples 

For each halo, we create a tracer of the DM consisting of 
10 6 randomly sampled DM particles. To constrain the po¬ 
tential profile of a halo all the way out to the virial radius, we 
adopt an outer cut of 300 kpc, which is slightly larger than 
the virial radius of the Aquarius haloes (200 to 250 kpc). 
We also adopt a n inner cut of 1 kpc , chosen to avoid con¬ 
vergence issues ([Navarro et al.l 1 20 id ). and to suppress the 
effect of any ambiguity in the definition of the centre of a 
real halo. By default, we use all the particles within the 
above radial range, no matter whether the particle belongs 
to the Friends-of-Friends halo or not. The Hubble flow is 
ignored throughout this analysis, since the scale at which it 
becomes important is given by GMJR ~ (HR) 2 /2, yielding 
R ~ 1 Mpc /h for a MW sized halou 

1 We have explicitly checked that including the Hubble flow pro¬ 
duces little difference in our results. 


3.2 Tagged Star Samples 

In reality one does not, of course, observe dark matter di¬ 
rectly. A realistic tracer population would be the stars in the 
halo of a galaxy. In this section we apply the oPD F method 
to th e Aquarius stellar haloes calculated by ICooper et al.l 
(l 201 (tl . These stars are identified in the output of the dark 
matter only simulation by tagging dark matter particles 
over time following the star formation history given by 
the GALFORM semi-analytical model of g alaxy forma¬ 
tion ([Cole et al.lll994l . r200d : lBower et al.ll2006h . The dynam¬ 
ics of the stars are then identical to the dynamics of the 
tagged dark matter particles. Since the dark matter par¬ 
ticles are dissipationless, this tagging method does not re¬ 
solve stellar discs. Nor does it take into account the effects 
of baryon dissipation on the gravity of the system. As a 
result, the distribution of stars in the inner galaxy is not 
quite realistic. Despite this limitation, the particle-tagging 
method provides a realistic model for the stripping and dis¬ 
tribution of accreted stars in the simulated outer halo, since 
the accreted stars follow the same collisionless d ynamics as 
the d ark matter particles on large scales fsee lLe Bret et all 
120151. for a controlled comparison of p article-tagging to hy - 
drodynamical simulations). Recently, ICooper et all ( 2013 ) 
have applied this technique to large-scale cosmological sim¬ 
ulations and have shown that it produces galactic surface 
brightness profiles that agree well with the outer regions of 
stacked galaxy profiles from SDSS. 

To test the oPDF method with a realistic tracer pop¬ 
ulation, for each halo we use the accreted stars from the 
particle-tagging technique. In addition, we exclude particles 
inside 10 kpc of each halo as the presence of a disc in a 
real galaxy violates the spherical symmetry assumption for 
the potential, and because the lack of such a disc in the 
simulated halo makes the mock data less realistic at small 
radii. As with the dark matter tracers, an outer radius cut 
of 300 kpc is applied to each halo. In a forthcoming paper 
(Wang et al., in prep), we will extend this study to a larger 
sample of Local Group haloes in which the stars are taken 
from hydrodynamical simulations. 

Due to the limited resolution of the dark matter simu¬ 
lation, each tagged particle may represent many stars with 
varied stellar masses, and one dark matter particle may be 
tagged multiple times representing stars formed at different 
epochs. However, the dark matter particles in the original 
simulation are followed dynamically without knowledge of 
the stellar mass weighting or multiple tagging. Hence the 
dynamics of these tracers are only resolved to the level of 
the tagged dark matter particles □ For the purpose of dy¬ 
namical modelling, we mainly use the unique set of tagged 
particles without any stellar mass weighting. This leaves us 
with 5 - 8 x 10 5 unique tagged particles for each halo in the 
level-2 simulations. In the following, we continue to use the 
term stars to refer to these unique sets of tagged particles. 


2 From a statistical point of view, the weighted distribution con¬ 
tributes an additional uncertainty to the stellar mass of each par¬ 
ticle, making the star particle counts in bins a Compound Poisson 
process rather than a Poisson process. So strictly speaking, the 
current likelihood model does not apply to the weighted distribu¬ 
tion. 
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Figure 1. The mass profile (scaled as M(< R)/R) of halo A. In 
the upper panel, the grey shaded line shows the true mass profile 
of the dark matter distribution, while the different coloured lines 
show NFW profiles from maximum likelihood fits within 50,100 
and 300 kpc respectively. The black dashed line (labelled virial) 
shows Eq. CD- It crosses each coloured line at the virial radius of 
each profile. The vertical dotted line marks the scale radius, r s . 
The lower panel shows the ratio of the fitted mass and the true 
mass as a function of the enclosing radius. 


3.3 Template profiles: defining the true potential 
and halo parameters 


To fit the halo potential using the oPDF method and assess 
any biases in the fit, we need to parametrize the potential 
with some functional form and also define the true parame¬ 
ters of the potential function. 

On e choice of parametrizatio n is the widely used NFW 
profile ([Navarro et al.|[l996Ul997h . 


p(r) = 


_ Ps _ 

(r/r s )(l + r/r s ) 2 ’ 


( 10 ) 


where r s is the scale radius at which din p/d In r = —2 and 
p s sets the density at this radius. These two parameters can 
be analytically related to the virial mass and concentration 
parameters. The virial mass is defined as the mass inside a 
virial radius, R v , where the enclosed density is A v times the 
critical density of the universe 

M = ^-A v p c lfi. (11) 

Throughout this paper, we adopt A v = 200. The con¬ 
centration parameter is defined as c = R v /r s . With this 
parametrization, one may choose to use the best-fitting pa¬ 
rameters of the density profile as the true parameters. How¬ 
ever, such a choice would be problematic if an NFW profile 
is not a good description of the halo density profile in ques¬ 
tion, in which case the best-fitting NFW parameters could 
depend on how the fit is performed. To demonstrate this, we 
fit the density profile of halo A using a maximum likelihood 


method. Note that dynamical modelling is not involved here, 
and the ht is purely to characterize the truejnass distribu¬ 
tion of the halo. The extended likelihood ([Barlow! Il990h . £, 
can be written as 

In C = ^2 In p(r,) ~ N pied , (12) 


where N pre d = J window p( r )/ m p ^ 3f is the predicted number 
of particles in the data window, with m p being the particle 
mass, and p(r) the NFW density profile given by Eq. (H0l) 
with parameters (p s ,r s ). ri is the radial coordinate of the 
z-th particle and the summation runs over ah the particles 
in the data window. This method is, in the limit of infinites¬ 
imal bins, equivalent to fitting to a binned profile provided 
one takes account of the Poisson distribution of the counts 
inside each bin. We ht the dark matter distribution around 
halo A over several different radial ranges, with outer cuts 
of 50,100 and 300 kpc respectively. The best-fitting mass 
profiles along with the real mass profile are shown in Fig. [I] 
It is obvious that the fits differ from each other, and none 
of them describes well the full mass profile out to the virial 
radius. The inferred virial masses can differ by more than 
30%. We note that halo A is an extreme example which de¬ 
viates grossly from NFW, while the remaining four Aquarius 
haloes agree much better with the NFW form. 

Given the poor performance of the NFW parametriza¬ 
tion for halo A, it would be problematic to define the true 
halo mass, concentration or potential parameters from a 
best-fitting NFW profile. Put another way, any ht that 
adopts an NFW parametrization also suffers from systemat- 
ics introduced by deviations of the real halo profile from the 
NFW form. To eliminate this systematic uncertainty, we will 
describe the potential using parametrized template profiles 
that are able fully to match the true profile. For each halo, 
we first extract the true potential profile from the spheri¬ 
cally averaged density profile. Specifically, the potential at 
a given point is evaluated as 

-^(r) = G^^ + G]r^, (13) 

ri<r r^r 

where n and rrii are the radial position and mass of the z-th 
particle. In practice, the profile is extracted at a sequence 
of radii and then interpolated at any other radius. Once a 
true profile is extracted, we generalize it to a two parameter 
family by varying its scale and amplitude. Specifically, for 
each real profile, z/>(r) = /(r), we generate a parametric 
template as 

ip(r) = Af (^ , (14) 

where A and B are dimensionless scale parameters. These 
two parameters can be mapped to M and c following the 
procedure in Appendix [Al The true parameters (Mo,co) of 
the halo are unambiguously defined by locating where in 
the true density profile the spherical overdensity matches 
the virial overdensity criterion and where the profile has a 
logarithmic slope of —2. 

We will consider both the NFW and the template 
parametrizations when fitting the potential. 
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Figure 2. The phase space distribution of particles in halo A. Top and bottom panels show the distributions in E — 6 and L — 6 spaces 
respectively. The left column shows the distributions of all the particles in the sample. The middle column shows that with subhalo 
particles removed. The right column shows the distributions of subhalo particles alone. Each subhalo with more than 1000 particles is 
also marked by a red circle in the right hand panels. In all the panels, only particles with 1 < r < 300 kpc are used. The white lines are 
the median 6s. The images colour code the number of particles in each pixel. The contrast of each panel has been individually optimized. 


4 APPLICATION TO DM HALOES 

4.1 The dynamical state of Aquarius haloes 

Once the real potential is known fEa. 1131) . we can examine 
the distribution of particles in (0, E, L) space prior to any 
fit. According to Eq. ©, for any system in a steady state, 6 
should be uniformly distributed for particles in any bin of E 
or L. In Fig. [2] we show the example of halo A in such coor¬ 
dinates. In the left panels, all the particles within 1 — 300 kpc 
from the halo centre are used. We are not concerned with 
the distributions along the E and L directions. Along the 
6 direction, overall, at fixed E or L the particle distribu¬ 
tions are close to uniform. However, one can still identify 
clumps in phase space which perturb the uniformity. In the 
rightmost panels on ly particles from subhaloes identified by 
subfind ([Soringel et aklEoP li are plotted. The coordinates 
of subhaloes with more than 1000 particles are overplot¬ 
ted as red circles, with larger circles corresponding to more 
massive subhaloes. The remaining particles representing a 


smooth component are plotted in the middle column. Com¬ 
paring the three columns, it is obvious that substructures 
introduce perturbations to the uniform 0-space distribution 
of the host halo. These perturbations are twofold: first the 
particles inside subhaloes are locally clustered and break the 
uniform distribution; secondly, the potential of the subhaloes 
exists as perturbations to the potential of the smooth host 
halo, affecting the orbits of nearby particles. The existence 
of locally clustered structures makes the real particle distri¬ 
bution noisier than a Poisson realization of a smooth uni¬ 
form field, and degrades the consistency between the two. In 
principle, substructures can be defined as locally overdense 
structures in phase space, and a phase space substructure 
finder could be designed to excise them and optimize the 
uniformity of the distribution of the remaining background 
particles. In practice, as we find subhaloes with subfind, 
the removal of substructures may not always increase the 
dynamical uniformity of the system unless the substructure 
finder is designed to do so. 
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Compared with the smooth component, subhaloes oc¬ 
cupy relatively low binding energy and high angular mo¬ 
mentum orbits. Despite the clumpy distribution of subhalo 
particles, they do not appear to bias the 0 distribution sig¬ 
nificantly in any particular direction. We will come back to 
this point when fitting the mass and concentration of the 
haloes. 

In Fig. O we explore deviations from a steady state of 
the DM tracer at different values of r, E and L in terms of 
the normalized mean phase deviation, 0, which measures 
the discrepancy level from a uniform distribution. For each 
halo, we calculate the mean phase within bins of phase-space 
coordinates r, F or L. We create the bins with equal numbers 
of particles per bin, so that they have the same statistical 
noise, allowing direct comparison of © across the bins. The 
bins are labelled by the percentiles of the respective sorted 
phasespace coordinate, r, F or L. Recall that, if the tracer 
is in a steady state, then in the large sample limit, 0 is 
distributed like a standard normal variable. 

Consistently with the physical picture displayed in 
Fig. El the DM particles have a mean phase deviation 
broadly consistent with zero. As seen from the left column, 
the discrepancy is most significant at large radius, low bind¬ 
ing energy and high angular momentum, revealing a higher 
level of systematics at these locations. Note low E and high 
L regions are also where subhaloes are most abundant as 
seen in Fig. [2] and it is also well k nown that subhaloes tend 
to occupy the outer halo (see, e.g.. lSpringel et ai1l2Q 08h The 
panels of the right hand column are the same as the left, but 
with subhalo particles removed from the tracer. In calculat¬ 
ing the radial profile, the radial limits of the data window, 
t min and r ma x, have been adjusted to the bin edges, so the 
radial profile examines the local uniformity of particles. Af¬ 
ter removing subhalo particles, local dynamical consistency 
is significantly improved at large radius. However, we see 
from the middle and bottom panels of Fig. [3] that this has 
little effect on the dynamical consistency within individual 
orbits over the full radial range. 

Note that 0 correlates with the depth of the proposed 
potential and a positive 0 indicates the current potential 
is deeper than a best-fitting potential (Paper I). As seen 
in Fig. [3j at large r, L and low E, the mean phase de¬ 
viation can be significantly higher than one would expect 
from a uniform distribution, which would lead to a level of 
systematic uncertainty significantly larger than the statisti¬ 
cal noise in the best-fitting potential. However, overall the 
fluctuation is still stochastic with no preferred sign. This 
indicates that if one is going to fit the potential, then devi¬ 
ations from our model assumptions are unlikely to bias the 
model parameters in a particular direction; instead, the bi¬ 
ases would fluctuate stochastically. Despite this, these biases 
are still systematic rather than statistical in nature, as they 
are tied to the model assumptions, not to the sample size. 
In the following section, we aim to quantify the level of such 
systematic uncertainty in the best-fitting parameters of the 
halo potential. 


4.2 Fitting the halo potential with DM as tracers 

With the potential functions and their true parameters de¬ 
fined in Section [3731 we can proceed to fit the potential pro¬ 
files with our oPDF method, and quantify the level of sys- 



1.00 1.02 1.04 1.06 1.08 1.10 
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Figure 4. The 1, 2 and 3cr confidence contours for the full sample 
of halo A fitted with the oPDF likelihood using the template 
profile. The parameters are in units of their true values. 


tematic uncertainty in the fitted parameters. We adopt the 
binned radial likelihood estimator, with 50 logarithmic bins. 
Fitting with a different number of radial bins gives consis¬ 
tent results@For both NFW and template parametrizations, 
we fit two datasets: 1) all the dark matter particles inside 
1 — 300 kpc, i.e, the full sample; 2) the former but with all the 
subhalo particles removed, i.e, the smooth sample. Because 
we aim to quantify the systematic uncertainties due to devi¬ 
ations from model assumptions, we need to make sure that 
the statistical noise, which is determined by sample size, is 
small enough. As an example, in Fig. [4] we show the sta¬ 
tistical confidence contours of halo A from the template fit. 
These error estimates are consistent with the scatter among 
independent subsamples of the parent halo. The la error is 
around 0.005 dex for our sample of 10 6 particles, quite con¬ 
sistent with our expected scaling of 0.1/a/ iV/1000 dex. Such 
an accuracy should be sufficient for detection of systematic 
biases larger than 1%. 

The best-fitting parameters in units of the true parame¬ 
ters are plotted in the left panel of Fig. [5] Overall, the fitted 
(M, c) parameters largely agree with their true values, with 
a bias generally smaller than 10%. The typical bias quanti¬ 
fied by the scatter among the five haloes is ~ 5% as listed in 
Table [3] For each parametrization and dataset, we combine 
the five haloes to estimate a mean and a covariance matrix 
for the parameters, and plot the one-sigma contour for a bi¬ 
variate Gaussian with the estimated mean and covariance. 
Note these contours are an estimate of the systematic uncer¬ 
tainties, since the statistical noise of the model is negligible 
given the sample sizes. Consistently with our expectation 
from the mean phase profiles, there is not a definitive sys¬ 
tematic bias but rather, as far as we can tell from the small 


3 Adopting t he Anderson-Darling estim ator described in Pa¬ 
per I (see also iBeloborodov &; Levinll20043 increases the parame¬ 
ter scatter to ~ 20%, due to its poorer accuracy. 
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Figure 3. The normalized mean phase deviation profile of Aquarius DM haloes. From top to bottom, we bin the halo particles according 
to their r, E , L coordinates respectively, with equal numbers of particles in each bin. The mean phase deviation, ©, is evaluated in each 
bin, and plotted as a function of the percentile values in the respective coordinate. Different colour lines represent different haloes. The 
dashed and dotted reference lines mark the 0 and =b3cr discrepancy levels. The left panels show the profiles of the full sample. The right 
panels show that with subhalo particles removed. See the online version for a coloured plot. 


sample of haloes, the scatter is mostly stochastic from halo 
to halo. 

Comparing the fits with and without subhalo particles, 
there is not a significant improvement in the latter. When 
NFW profiles are adopted in the fits, the accuracy is compa¬ 
rable to that achieved with template profiles in most cases. 
This reflects the fact that most haloes are well described by 
NFW profiles. Note that the significantly larger confidence 
regions in NFW fits as marked by the ellipses in Fig. [5] is 
caused purely by halo A, whose dynamical fit shows a bias 
in concentration up to 30%. This is due to the fact that 
the density profile of halo A differs significantly from NFW, 
as is evident from Fig. \l\ In Fig. [6] we show the disagree¬ 
ment in halo A from a different perspective, by comparing 
the NFW parametrization of the halo potential with the 
true potential. When the set of true halo parameters are 
used, the NFW potential is consistently overestimated in¬ 
side the halo. The dynamical fit adjusts the parameters so 
that the NFW potential agrees with the true potential to 
within 5 percent for most of the radial range. The best-fit 


NFW potential agrees better with the true potential, reflect¬ 
ing the fact that the dynamical fit largely recovers the true 
potential by force-fitting the NFW parametrization, despite 
giving different parameters from the true values. It is quite 
interesting to see that when the template profile is adopted, 
halo A does not appear to be more biased than the other 
haloes, meaning that a deviation from the NFW form does 
not necessarily mean a lack of equilibrium. 

In the right panel of Fig. [5j we compare our fits to 
those obtained from a conventional DF method that de¬ 
scribes the phase space density onl y as a function of (E, L ) 
of the particles (|Wang et al.ll201a ). Specifically, the phase 
space probability is assumed to have the form d P(f,v) = 
f(E)L~ 2(3 d 3 r d 3 v, where [3 is a parameter describing the ve¬ 
locity anisotropy, with f(E) further determined by inverting 
a double-powerlaw tr acer density profile inside an NFW po¬ 
tential (see Eq. 12 in I Wang et alJEoi5l for further details). 
This distribution function describes a family of models with 
constant anisotropies, while in gen eral more flexible mod- 
els can be constructed (e.g. IWoitak et ahll2008l : IPosti et al.l 
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Figure 5. Left: fitted parameters of Aquarius haloes using the radial likelihood estimator. Different shaped symbols denote different 
haloes. The red and blue colours denote the fitting results using NFW and template profiles (labelled TMP) respectively. In both cases, 
the open symbols show the fits with the full samples, while the filled ones show those for the smooth samples, i.e., with subhalo particles 
excluded. For each combination of sample and profile, we combine the five haloes to estimate a mean and a covariance matrix for 
the parameters, and plot the la contour (the ellipses, open or filled) in the same style fo r a bi variate Gaussian with the estimated 
mean and covariance. Right: same as the left, but also showing the fits from I Wang et al.1 (l2015h to the smooth DM sample using a 
f(E,L) = L~ 2 PF(E) model (green symbols and ellipse). See the online version for a coloured plot. 



Figure 6. The potential profile of halo A and its NFW 
parametrization. We plot the ratio between the parametrized 
NFW potential, ^» and the true potential, t/j o, of the halo as 
a function of radius. For the dynamical modelling only the po¬ 
tential difference is relevant and so the zero-points of the NFW 
potentials have been adjusted to produce the true potential value 
at the virial radius. The red solid line corresponds to the NFW 
potential profile using the true parameters and the green dashed 
line to the NFW parameters found from the oPDF likelihood es¬ 
timated from all the dark matter particles. 


120151 : [Williams &; Evansll2015bh . In contrast to the fairly un¬ 
biased fits with our method, this f(E : L) method suffers 
from a ^ 50 % net bias in the parameters. As discussed in 
IWang et al.l (l2015h . this can be attributed to the fact that 
the f(E , L ) DF only describes gravitationally bound sys¬ 
tems by construction (see also section 6.3.1 of Paper I), and 
struggles to match the distribution of the loosely bound par¬ 
ticles. Because our oPDF method has no prerequisite on the 
distribution of orbits (hence no prerequisite on the energy 
distribution), our fits show no such net bias. On the other 
hand, the f(E , L) fits exhibit a comparable amount of halo 
to halo scatter in the parameters to ours, reflecting that our 
method does capture the minimum irreducible uncertainty 
associated with steady-state models. 


5 APPLICATION TO MOCK STELLAR 
HALOES 

In Fig. 0 we show the phase space distribution of stars in 
halo A. Unlike the DM tracer which is only slightly per¬ 
turbed by subhaloes, the halo stars are dominated by those 
in the satellite subhaloes. The mass fraction contained in 
satellite galaxies is 50% — 70% (Table |T|) in the radial range 
of interest. These satellite stars are obviously not in equi¬ 
librium with the rest, and can be observationally identified 
and removed as satellite galaxies. In what follows, we will 
only use the “smooth” component of halo stars, i.e, those 
excluding satellite stars, as our tracer sample. In total, each 
halo has (2 — 5) x 10 5 smooth star particles within 10 — 300 
kpc, yielding a statistical uncertainty of ~ 2% in mass and 
concentration. 

The mean phase deviation profile is shown in Fig. [8] As 
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Table 1. Basic properties of the stellar haloes, within 10 — 300 kpc. Ntot is the total number of star particles, N srnt h/Ntot is the fraction 
of particles in the smooth component and m smt h/'/ntot is the mass fraction of particles in the smooth component. 


Halo 

Wot 

N sm th/Ntot 

^smth/^tot 

D 

A 

5.4 x 10 5 

0.50 

0.34 

9.6 

B 

7.4 x 10 5 

0.72 

0.36 

9.6 

C 

5.3 x 10 5 

0.70 

0.43 

5.9 

D 

8.1 x 10 5 

0.60 

0.52 

5.5 

E 

5.1 x 10 5 

0.47 

0.18 

9.5 


for the profile of DM tracers, overall 0 is consistent with 
zero, with the highest scatter seen at large radius, low bind¬ 
ing energy and high angular momentum. Note that stars 
with low bindin g energies are also those that have been ac¬ 
creted recently (|Wang et al.ll2015h . The scatter in the star 
profiles also appears higher than that in the DM case. 

The fits using stars are plotted in Fig. [9] for both the 
template and NFW profile models. For individual haloes, the 
deviation from the true parameters can be as high as 40%. 
For comp arison, the fits and la contours from the f(E, L) 
method of I Wang et al.l (l2015h and those from the DM tracers 
in the previous section are also plotted. 

Overall, we do not observe a statistically significant 
net bias in the fits with the current sample of five haloes, 
even though the f(E , L ) method applied to stars is only 
marginally unbiased at the lcr level. In other words, the sys¬ 
tematic bias varies from halo to halo in a stochastic way. 
Despite this stochastic behaviour, we have checked that the 
systematic bias does not change with sample size, so it is 
indeed a systematic rather than a statistical error. There is 
a negative correlation between the mass and concentration 
parameters in the template fits, which is similar to the cor¬ 
relation in the statistical noise of the two parameters. This 
correlation is absent in the NFW fits only because halo A 
is not well described by NFW and this biases the fit signifi¬ 
cantly. A viable explanation of this behaviour of the system¬ 
atic bias lies in the deviation of the tracer population from a 
steady state. For example, the existence of correlated phase 
angles in streams and caustics implies that different tracer 
particles are not independent. As a result, the constraining 
power of a set of particles in a stream is less than that of an 
equal number of independent particles. In the large sample 
limit, when each stream is sufficiently sampled, the errors on 
the inferred model parameters do not vanish but are limited 
by the effective number of independent streams or particle 
clumps. This is an intrinsic property of each halo. Hence it 
is understandable that we are left with irreducible stochas¬ 
tic biases in well sampled haloes. In addition, these residual 
errors are expected to exhibit similar parameter correlations 
to the statistical noise. Note that while the DM fits exhibit 
only ~ 5% scatter, the scatter for the stellar fits is typically 
- 20%. 

Since stars and DM tracers have different E-L distribu¬ 
tions, they occupy a statistically different set of orbits. This 
implies the tracers potentially have different spatial distri¬ 
butions. As a result, they could be sampling different parts 
of the halo, or the same region but with different weights 



Figure 8. The mean phase deviation profile of Aquarius stellar 
haloes. This is the same as the right hand side of Fig. [3] but for 
star particles. From top to bottom, we bin the star particles ac¬ 
cording to their r, E and L coordinates respectively, with equal 
numbers of particles in each bin. The mean phase deviation, ©, is 
evaluated in each bin, and plotted as a function of the percentile 
values in the respective coordinate. Different coloured lines rep¬ 
resent different haloes as indicated in the legend. The dashed and 
dotted reference lines mark the 0 and d=3cr discrepancy levels. 
Only the smooth component of the halo stars is used. See the 
online version for a coloured plot. 
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Figure 7. As Fig. [2] but for the stars in the level 2 halo A. From left to right the distributions of all the stars, the “smooth” component 
(all stars excluding satellites), and those in satellite subhaloes. The contrast of each panel is individually optimized. 


given to the local deviations of the halo potential. It is pos¬ 
sible that the different sampling has resulted in the stars 
yielding a large scatter in the inferred halo properties. To 
see whether this is the case, we select dark matter parti¬ 
cles that have the same E-L distribution as the stars to 
create a star-like dark matter sample. Subhalo particles are 
removed from the DM and star samples before sampling 
the E-L distribution and the radial coordinates are ignored 
when constructing the samples. The same fitting procedure 
is then applied to this star-like dark matter sample. The re¬ 
sults are shown in Fig. [10] By drawing a sample with the 
stellar P(E, L) from DM particles, the scatter in the fits 
is actually slightly decreased (probably due to the removal 
of the less virialized outer halo particles) compared to the 
DM fits, and is much smaller than that in the star fits. This 
shows that the tagged stars are indeed in less of a steady 
state than the dark matter. 

We remark that even though for the star samples the 
scatter in Fig. [8] is higher at large radii, low binding energies 
and high angular momentum, we do not detect a systematic 
decrease in the biases of the fitted parameters as we exclude 
the large radii, low energy or low angular momentum re¬ 


gions. The lack of systematic improvements in bias when 
excluding the regions with large scatter is consistent with 
our previous argument that the biases are limited by the 
effective number of independent streams. Also note the bias 
and scatter are separate quantities, and we have observed 
that in the high scatter regions the mean profile does not 
appear more biased. 


6 HALF MASS CONSTRAINT 

In Paper I we used Monte-Carlo samples to demonstrate 
that the mass profiles are best constrained near the median 
radius of the tracer population. Similar best-constrained 
masses also exist in se veral previous stud i es using very dif¬ 
ferent methodologies ([Walker et ahl 120091 : IWolf et all 1 20 id : 
lAmorisco &; Evanal2011 . see section 5.1 of Paper I for more 
discussion). Here we revisit this discovery with the Aquar¬ 
ius haloes. In Fig. [Tl] we plot the constrained mass profile 
from the DM and star tracers in halo A. Note our likeli¬ 
hood method constrains not only the characteristic mass, 
but also the shape of the mass profile. If the profile shape 
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Figure 9. As Fig. [5] but showing the results of the dynamical fits 
to the halo stars. As indicated in the legends, symbols of different 
shapes represent different haloes, while different colours distin¬ 
guish different datasets and model profiles. Red and green show 
fits to stars adopting NFW and te mplat e profiles respectively; 
blue shows the fit from I Wang et all (l2015h to stars combining ra¬ 
dial and tangential velocities using a specific f(E , L) model; grey 
shows the template fit applied to the smooth DM tracer. The 
symbols are the results of fits to individual haloes, while the large 
ellipses mark the estimated 1 -cr confidence regions for each type 
(i.e., combination of dataset and model) of fit estimated from the 
sample of five haloes. See the online version for a coloured plot. 



0.6 0.8 1.0 1.2 1.4 

m/m q 


Figure 10. Comparison between the results of stellar and the 
DM dynamical fits. The “Star” and “DM” fits are the same as in 
Fig. [9] using stars and DM particles respectively. The “dmStar” 
fit uses a sample of DM particles selected to have the same E- 
L distribution as that of the stars. Subhalo particles have been 
removed in all three cases and only the template profiles are used. 
See the online version for a coloured plot. 


Figure 11. The template-fitted mass profile of halo A, using 
smooth DM (red) and smooth star (grey) particles as tracers re¬ 
spectively. The shaded regions are the 1-cr constraints on the mass 
profile, normalized by the true profile. The vertical lines mark the 
half-mass radius of the two tracers. The green dotted line is the 
best-fitting profile from the star-like DM tracer (i.e., “dmStar” in 
Fig. [TO]). 


is biased then the bias in enclosed mass varies with radius. 
Consistent with our previous findings, the mass is best con¬ 
strained near the half-mass radius of the tracer. Comparing 
these best-constrained masses, the bias in stars is still signif¬ 
icantly larger than that in DM. This is also consistent with 
our test using a star-like DM tracer in Fig. 1101 where we find 
that the different samplings of the star and DM tracers are 
not the cause of the different bias levels. For comparison, the 
best-fitting profile of the star-like DM tracer is also plotted, 
and shows a bias comparable to the original DM fit. As listed 
in Table [2] stars yield an average bias of ~ 5% at r i/ 2 when 
using the template fits, while the DM yields only ~ 1%. The 
star-like DM tracers have 7 * 1/2 close to that of the stars, but 
gives almost no bias at 7*i/ 2 . As far as the constraints at r*i/ 2 
are concerned, fitting with NFW profiles gives quite similar 
results to template fits, indicating that the half mass con¬ 
straint is less sensitive to the adopted functional form for 
the halo density profile (see also Paper I). However, models 
with extra assumptions could still lead to significant bias at 
7 * 1 / 2 • For example, fitt i ng th e DM tracers with the f(E,L) 
method in I Wang et al.l (l2015h produces an average mass bias 
of 13% at 7 *i/ 2 (Table [3j. 

We emphasize that since we are not only interested in 
the mass constraint at a single radius but also in the full pro¬ 
file, any parametrization of a specific density profile should 
be equivalent. As long as the constraints are fully described 
in terms of the parameter covariance or the 2-dimensional 
confidence contour, the constraints on the full profile can al¬ 
ways be recovered and translated to constraints in any other 
parametrization of the profile. Our parametrization is inten¬ 
tionally chosen to constrain the most popular parameters, 
the virial mass and concentration of haloes. 
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Table 2. Tracer half-mass radius 7*1/2 an d mass bias b^/2 at 7*1/2, from fits to the smooth component of DM, star, and “dmStar“ tracers. 
“dmStar” refers to DM tracers selected to have the same E-L distribution as the stars, as in Fig. 1101 We list the biases from template 
fits by default. For DM and star tracers, biases from NFW fits are also given in parenthesis. 


Halo 

DM r*!/ 2 /kpc 

DM bij 2 

Star 7^/2/kpc 

Star bi /2 

dmStar r 1 f 2 

dmStar b 1 j 2 

A2 

103 

0.00 (0.05) 

41.7 

-0.07 (-0.02) 

41.6 

-0.01 

B2 

85.4 

0.01 (0.03) 

18.8 

0.05 (0.05) 

19.0 

0.01 

C2 

86.2 

0.03 (0.05) 

48.2 

0.03 (0.04) 

47.4 

0.00 

D2 

103.1 

-0.02 (0.00) 

32.8 

0.04 (0.05) 

33.1 

0.00 

E2 

90.1 

-0.00 (0.00) 

18.6 

0.04 (0.02) 

18.6 

0.00 


Table 3. Summary of the different fits to the halo density profile. For each combination of data and method, we list the fitted parameters 
averaged over the five haloes (x) and their halo-to-halo standard deviation (cr) in the form x d= g. The mass (M) and concentration (c) 
parameters are normalized by their true values, Mo and co- The mass bias at the tracer half-mass radius, &i/ 2 , is also listed in the 
same form. Different columns refer to different combinations of data and methods. “DM-Full”, “DM” and “Star” refer to full DM, the 
smooth DM (DM-Full excluding subhalo particles), and smooth star tracers. “dmStar” refers to DM samples selected to have the same 
E -L distribution as s tars. “NFW” and “TMP” refer to fits using NFW or template potential profiles. f(E,L) refers to the (r,v r ,vt) fit 
in I Wang et alii i2015l i using an f(E , L) distribution function. 



DM:NFW 

DM:TMP 

DM-Full:TMP 

DM :f(E,L) 

StariNFW 

Star:TMP 

Star:/(£, L) 

dmStar:TMP 

M/M 0 

0.97 ±0.07 

1.00 ±0.05 

1.02 ±0.04 

1.54 ±0.05 

1.00 ±0.25 

1.10 ±0.23 

0.82 ±0.16 

0.99 ±0.03 

c/c 0 

0.99 ±0.17 

0.98 ±0.06 

0.96 ±0.04 

0.48 ±0.13 

0.99 ±0.16 

0.97 ±0.18 

1.12 ±0.14 

1.01 ±0.03 

&1/2 

0.02 ±0.02 

0.00 ±0.01 

-0.01 ±0.02 

0.13 ±0.02 

0.03 ± 0.02 

0.02 ± 0.04 

-0.02 ±0.01 

0.00 ±0.01 


7 SUMMARY AND CONCLUSIONS 

We have applied our oPDF estimator to tracers in five simu¬ 
lated haloes from the Aquarius project, to study the level of 
systematic biases in fitting the mass distribution of Milky- 
Way sized haloes. We focus on effects from the parametriza- 
tion of the halo potential, the existence of subhaloes, and the 
types of tracer used. Assuming a spherical symmetric poten¬ 
tial, our method only makes use of the dynamical equilib¬ 
rium of the tracer. As a result, the level of systematic biases 
detected in our analysis can, in general, be interpreted as 
the minimum level of bias present in any time-independent 
DF modelling of dynamical tracers that assumes a spherical 
potential. With our sample of five haloes, we do not have a 
reliable detection of a common bias in our method towards 
any particular direction in parameter space. Instead, we fo¬ 
cus on characterizing the average amplitude of the bias in 
each fit. We quantify this as the rms scatter of the biases 
for individual haloes, and summarize them in Table El The 
method works very well on DM tracers, with a level of sys¬ 
tematic bias at only ~ 5%. Assuming an NFW profile does 
not significantly affect the fits in most cases, except in one 
case out of five where the density profile of the halo (Aq-A) 
differs significantly from the NFW form, leading to a much 
larger bias (~ 30%) when adopting the NFW profile. How¬ 
ever, the deviation from the NFW profile in halo A does not 
affect the equilibrium of the DM tracer. Subhaloes exist as 
perturbations that give rise to deviations from steady state 
for the tracers, but only affect the dynamical fits using DM 
tracers by ~ 1%. In contrast to the fairly good fits to the 


DM tracers with our metho d, a conventional DF fit adopt¬ 
ing a specific f(E,L) DF (|Wang et al.1 12015 ) yields a net 
bias of 50% in mass and concentration on average. This is 
caused by the additional assumptions made in the f(E,L) 
DF that restricts the allowed distribution of orbits. On the 
other hand, the halo-to-halo variation in the bias is compa¬ 
rable to that in our method, demonstrating that our method 
gives the minimum uncertainty in mass modelling assuming 
time-independent DFs. 

Applying our method to mock stars results in a higher 
level of bias, b ~ 2 0%, comparable to that in the f(E , L ) DF 
method tested in I Wang et al.l 1201 5) which, however, also 
suffers from a non-zero net bias. The larger bias using star 
tracers is not due to different phase-space sampling by stars 
compared to DM particles: DM tracers constrained to sam¬ 
ple the phase space E-L distribution in the same way as the 
stars yield biases at the same level as the original DM trac¬ 
ers. The larger deviation of tagged stars from a steady state 
is not surprising because they involve the 1% most-bound 
particles of their host subhalo at the time of star formation. 
By definition, these particles are the most resistant to tidal 
stripping and subsequent mixing. Even though we only use 
the stripped population of tagged particles, they are still 
farther from equilibrium than the smooth component of the 
host halo. 

It is well known that dynamical tracers best- 
constrain the host mass n e ar th e tracer’s ha l f-mass 
radius, r-\ / 2 (IWalker et al.l 120091 : IWolf et ahl l20ld : 
lAmorisco &; Evansl 1201 lh ~ Although we adopt a vastly 
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different method from those analysis, a similar behaviour is 
also observed in our analysis. Near 7q/ 2 , the mass biases, 
61/2, are much reduced, and also become less sensitive to the 
functional form of the halo profile assumed in our model. 
Larger biases are still observed for stars, b 1 / 2 ~ 5%, com¬ 
pared with b 1 / 2 ~ 2% for DM tracers. The f(E,L) method 
that poorly fits the DM tracer produces a much larger 
bias, bi /2 ~ 10%. Although the bias at t*i / 2 is significantly 
smaller than the bias for the total mass, in reality bi/ 2 
together with the constraint on the profile shape at r*i/ 2 is 
equivalent to the joint constraint in the mass-concentration 
space. Given the full mass-concent rat ion covariance matrix, 
one can readily obtain the mass constraint at any radius 
including ri/ 2 . While ri/ 2 depends on the tracer, M and c 
are intrinsic properties of the halo. 

In t his work we freque ntly compare our results with 
those of IWang et al.l 1201 51 who used a specific model of 
the f(EjL) family to study the same haloes. We demon¬ 
strate that the extra assumptions in that model beyond 
time-independence and spherical symmetry have resulted 
in a worse performance compared with the oPDF. There 
are more flexible distri bution functions t hat can improve 
over the one assumed in I Wang et al.l (l2015ll. for example, by 
allowing for varying anis otropies (e.g. IWoitak et al.l [20081 : 
IWilliams fe Evansll2015bh . More generally, there may exist 
a true model that describes well the distribution function of 
the tracers. However, such a true model has to be known 
a priori to fit the tracers correctly, which is a highly chal¬ 
lenging task if at all possible. At the same time, any specif¬ 
ically proposed distribution function has generally has lim¬ 
itations stemming from extra assumptions over and above 
the Jeans theorem. These extra assumptions may not be 
obeyed by an arbitrary tracer sample. As such, the results 
obtained using the oPDF method which makes minimal as¬ 
sumptions are partic ularly robust , and t he comparison with 
the specific model of I Wang et al ] d201fih therefore serves to 
illustrate the limitations of restricted models. In particular, 
since the statistical noise has been controlled to be negligi¬ 
ble in our analysis, the level of systematic bias detected with 
the oPDF is the minimum level of systematic bias expected 
from any model that assumes a spherically symmetric and 
time-independent distribution function. 

Note that the stars used in this work are generated from 
a particle tagging method that is a relatively simple way of 
approximating the phase-space distribution of stars. Several 
factors, including insufficient mass resolution (the weighting 
and multiple tagging of star particles), the time discrete¬ 
ness of the tagging (the method works with snapshots), and 
the lack of dissipation and back-reaction on the potential 
from stars, could all potentially affect the degree of realism 
with which the tagged stars represent the dynamics of real 
stars. Due to these limitations, the results from the tagged 
stars should only be taken as indicative. Also note that only 
five haloes are studied in this work and these may not be 
very representative of our Milky Way halo. In a follow up 
work, we will apply the method to a larger sample of haloes 
modelled using SPH simulations of the Local Group for a 
more realistic assessment of the dynamical state of tracers 
in Galactic haloes. 
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APPENDIX A: PARAMETERS OF TEMPLATE 
PROFILES 

To connect template parameters A and B in Eq. (1141) to 
physical parameters, we can define 

A = '0s/'0 sO , (Al) 

B = r s / r s0 , (A2) 

where r s is a scale radius at which the profile has some prede¬ 
fined shape, and is the potential at r = 0. We choose r s to 
be the radius where din p/d In r = —2, to be consistent with 
an NFW parametrization. r s o and ^so are the corresponding 
quantities of the true profile. Hence this template profile is 
parametrized by (A, B) or equivalently by We can 

also define equivalent mass and concentration parameters. 
For each profile, the virial mass, M, and virial radius, R v , 
can be defined following the same spherical-overdensity def¬ 
inition as in Eq. CP, and the concentration can be defined 
through c = R v /r s , consistently with NFW0 The mass and 
concentration parameters of the true profile (i.e., the tem¬ 
plate with A = 1, B = 1), Mo and Co, are by definition the 
true parameters of the halo, and can be obtained unambigu¬ 
ously from the true profile without fitting. If the halo is per¬ 
fectly NFW, then the true parameters defined this way are 


also the best-fitting NFW parameters to the density profile. 
When the density profile differs from NFW form, however, 
the true parameters, Mo and Co, should be interpreted as the 
spherical overdensity mass and the contrast of the spherical 
overdensity radius, R v , to the slope —2 radius, r s , rather 
than being any best-fitting NFW parameters. 

With the template parametrization, the inversion from 
any set of (M, c) parameters back to (A, B) is also straight¬ 
forward. Note that the mass profile of the template scales 
as 

MM - (A3) 

where M(r) is the mass profile of the template with parame¬ 
ters (A, B ), ip f (r) is the derivative of the template potential, 
and m{r) is the true mass profile. Hence M — ABm(R v /B). 
After obtaining (R v , r s ) from (M, c), one can solve (A, B) as 
follows 


B = — 

TsO 


A = 


M 


Bm(R v /B) ‘ 


(A4) 

(A5) 

(A6) 


To create the templates numerically, we extract both the 
potential profiles and the cumulative density profiles p(< 
r) oc , 0 / (r)/r from the particle distribution of each halo. 
The p(< r) is provided to avoid the need for numerical dif¬ 
ferentiation of the potential profile. 


4 Although we have chosen r s to be the slope —2 radius, in prin¬ 
ciple r s can be defined to be the radius at any characteristic slope, 
with the concentration parameter being interpreted as the ratio 
between R v and r s . As long as the definition is consistent within 
the same template, the B parameter does not depend on the spe¬ 
cific definition of r s . 




